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Recent experiments on the Ba 3 XSb 2 09 family have revealed materials that potentially realise 
spin- and spin-orbital liquid physics. However, the lattice structure of these materials is complicated 
due to the presence of charged X^^-Sb^^ dumbbells, with two possible orientations. To model the 
lattice structure, we consider a frustrated model of charged dumbbells on the triangular lattice, with 
long-range Coulomb interactions. We study this model using Monte Carlo simulation, and find a 
freezing temperature, Tfrz, at which the simulated structure factor matches well to low-temperature 
x-ray diffraction data for Ba 3 CuSb 209 . At T = Tfrz we find a complicated “branching” structure of 
super exchange-linked X^^ clusters, which form a fractal pattern with fractal dimension df = 1.90. 

We show that this gives a natural explanation for the presence of orphan spins. Finally we provide 
a plausible mechanism by which such dumbbell disorder can promote a spin-orbital resonant state 
with delocalised orphan spins. 

PACS numbers: 75.10.Kt, 75.25.Dk, 75.47.Lx 


Recently there has been an intense search for materi¬ 
als exhibiting spin-liquid behaviour - materials beyond 
the “standard model” of condensed matter physics^. A 
particularly intriguing idea is of a spin-orbital liquid, in 
which not only the spin but also the orbital degrees of 
freedom remain fluctuating down to low temperature^“^. 

The Ba3XSb209 family, with 

.., has been shown 
to be a promising class of materials to realise spin-liquid 
behaviour. Ba3CuSb209 has been particularly well 
studied, and it has been suggested that the spin and 
orbital degrees of freedom associated with the Cu^+ 
ions form a spin-orbital liquid state^“^^. In the case of 
Ba3NiSb209, the pressure-synthesised 6 H-B structure 
has been proposed as an example of a spin -1 spin-liquid 
state^^”^^. 

An important starting point when trying to under¬ 
stand spin-liquid behaviour is knowledge of the lattice 
structure. In Ba3CuSb209 it has been suggested that 
the Cu^+ ions form a short-range honeycomb lattice^, 
and theoretical approaches have therefore concentrated 
on Cu^+ plaquettes formed of several hexagons^^’^^’^^. 
On the other hand, in the 6 H-B phase of Ba3NiSb209 it 
has been suggested that the Ni^+ ions form a triangular 
lattice^^. 

Here we argue that in neither case is this a good 
starting point for theoretical investigation, and instead 
one should consider a disordered “branch” lattice [see 
Fig. lb]. The evidence we present focuses in particular 
on Ba3CuSb2 09, but should be applicable to other mem¬ 
bers of the Ba3XSb2 09 family. Furthermore, we suggest 
that this type of correlated lattice disorder can promote 
spin-orbital liquid behaviour. 

In order to investigate the lattice structure of these ma¬ 
terials, we solve a frustrated model of interacting X^+- 
Sb^+ charged dumbbells [see Fig. 1 ]. We argue this is 
relevant to stoichiometric X=Cu, X=Ni in the pressure 
synthesised 6 H-B phase and potentially to pressure syn¬ 



FIG. 1: Charged dumbbells on the triangular lattice, (a) X^+- 
Sb^+ dumbbells of length z form a triangular lattice bilayer. There 
is an Ising degree of freedom associated to whether the dumbbell 
is orientated with X above Sb or vice versa. The equilibrium dis¬ 
tribution of dumbbells can be mapped onto a charge model, Tcoul 
[Eq. 1], which at low temperature orders in a stripe ground state 
[shown here], (b) Material realisations of Tcoul [Eq. 1] fall out of 
equilibrium at T = Tfrz, and the lattice structure can be studied by 
making simulations at this temperature. A snapshot of a typical 
lattice structure for X=Cu is shown, with blue and orange sites 
denoting different dumbbell orientations. Superexchange interac¬ 
tions link Cu^+ ions on dumbbells with the same orientation, and 
superexchange linked clusters are shown by blue and orange bonds. 

thesised X=Mn and X=Co. 

The X^+-Sb^+ dumbbells are surrounded by bioc- 
tahedra, and their constituent ions sit on the vertices of 
stacked triangular lattice bilayers^, as shown in Fig. la. 
Each dumbbell has two possible orientations with either 
the X 2 + or Sb 5 + on top. Electrostatically, the primary 
influence on the orientation of the dumbbells is the ori¬ 
entation of the other dumbbells - that is to say that the 
Ba^+, and remaining Sb^+ ions are electrostatically 
ambivalent as to the dumbbell orientation. 

This leads us to consider a Coulombic charge model, 
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where = ±1 is a normalised charge, i and j run over 
the sites of a bilayer triangular lattice [shown in Fig. 1], 
Vij = |r^ ~ i*j| and the charge distribution is constrained 
to have one positive and one negative charge on each 
dumbbell. We ignore interaction between dumbbells in 
different bilayers, and we provide a justification for this 
approximation below. 

In order to relate the charge distribution following from 
F^Coui [Eq. 1] to the lattice structure of the materials, it 
is necessary to understand the synthesis process. This 
is typically performed at high temperature (>1000°C), 
and the crystals are then slowly cooled to room temper¬ 
ature and below^’^. A characteristic timescale tcooi can 
be ascribed to this cooling process, and this should be 
compared to tfnp, the characteristic time for dumbbells 
to reverse their orientation. Close to the synthesis tem¬ 
perature, we assume that tfnp ^ tcooh and therefore the 
dumbbell orientation remains in thermal equilibrium as 
T is reduced. As the crystal is cooled, tfnp increases, and 
there is a temperature, Tfrz, below which tfnp ^ Aooi- la 
this regime the dumbbell dynamics is too slow to equilib- 
riate the system and the charge distribution is thus frozen 
in place. The dumbbell structure for any T < Tfrz can 
therefore be understood from studying the equilibrium 
dumbbell structure at T = Tfrz- 

The dumbbells in these materials are widely spaced, 
and one piece of evidence that they are dynamic at 
high temperature comes from the isostructural com¬ 
pound Ba 3 lrTi 209 ^^. Here the Ir-Ti dumbbells exhibit 
a markedly different low-temperature structure depend¬ 
ing on whether the material is slowly cooled from the 
synthesis temperature (1000°C) or quenched. 

This suggests a twofold strategy for understanding the 
lattice structure of these materials. 1) Simulate Tcoui 
[Eq. 1] as a function of temperature, and, by comparison 
with experimental data, determine the freezing temper¬ 
ature, Tfrz. 2) Simulate the model at Tfrz in order to 
extract detailed information about the lattice structure 
for all T < Tfrz- 

In order to simulate Tcoui [Eq. 1], it is first mapped 
onto an Ising model on the triangular lattice using Ewald 
summation^^. This leads to, 

-Ecoul ~ -Eq T — j ^ (2) 

id 

where = ±1 is an Ising spin, i runs over the sites of a 
triangular lattice and 'ipijiz) defines the interactions be¬ 
tween sites as a function of the dumbbell size, 2 ; [see Eig. 1 
and Appendix A]. Eor z ^ 0, Ecoui [Eq. 2] reduces to in¬ 
teracting Ising dipoles on the triangular lattice^^. Here 
we consider z = 0.46a as this is relevant to Ba 3 CuSb 2 09 ^. 

We have performed Monte Carlo simulations of Ecoui 
[Eq. 2] over a wide range of temperatures [see Ap¬ 
pendix B]. The ground state is 6-fold degenerate, and 
consists of alternating stripes of cr = 1 and a = — 1, 
parallel to either the A, B or C bonds [see Eig la].^^ At 
Tc/'ipnn ~ 0.19, with the nearest-neighbour interaction 
'0nn ~ 0.18 in the units of Eq. 1, there is an apparently 



FIG. 2: The dumbbell structure factor, S'(q), as predicted by sim¬ 
ulations of -Ecoul [Eq. 1]. The structure factor is plotted at a range 
of temperatures as a function of qx with qy = 27r/\/3 and qz = tt. 
From top to bottom: T/'ipnn = 0.45 (black), Tlip^n = 0.9 (blue), 
Tl'ipnn = 1.4 (green), Tj'ipnn = 2.4 (orange) and Tj'dnn = 3.2 (red). 
In the inset, the ratio R = S{0,27t/\^, qz)/S{27 v/3,27t/\^, qz) 
(small dotted arrow compared to large dashed arrow) is plotted 
as a function of T. 


1st order phase transition into a domain wall network 
state, as proposed in Ref. [30] for the Ising model with 
further neighbour exchange interactions. We postpone 
a detailed description of the low-temperature behaviour 
to another publication, and instead concentrate on the 
temperature region above the phase transition. 

For T > Tc we perform simulations to measure the 
dumbbell structure factor. In the absence of interaction 
between bilayers, this is given by. 


5'(q_L,gz) 
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where rx,i measures the position of dumbbell i in the 
plane of the triangular lattice. Here Qz = 27rzl/c, where 
for Ba 3 CuSb 209 the dumbbell height is z = 2.69A, the 
unit cell has a height c = 14.37A and I is measured 
relative to the structural Bragg peaks^. For / = 0 it 
is not possible to observe scattering from the dumbbell 
structure, as there is a destructive interference between 
X and Sb ions within the same dumbbell. Scattering is 
strongest when Qz = (2n + l)7r, where n is an integer, 
and for n = 0 this corresponds to I = cl{2z) ~ 3. Fig. 2 
shows S{qx, 27r/\/3^7r) at a range of temperatures, and 
there are diffuse peaks centred on = (=b27r/3, 27r/\/3) 
[see also Fig. 3a]. 

The diffuse nature of the peaks in S{q) [Eq. 3] corre¬ 
sponds to the absence of long-range order in the dumbbell 
structure. The width of the peaks at half maximum gives 
a measure of the correlation length, ^ig. Eor example, for 
T = 0.9V^nn [blue curve in Eig. 2], we find ^\s ^ 2a. In 
domains of this lengthscale the system is correlated in a 
stripe-like pattern [see Eig. 1]. 

The motivation for studying the dumbbell structure 
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FIG. 3: Comparison between the simulated structure factor fol¬ 
lowing from Fcoul [Eq. 1] and x-ray diffraction experiments for 
Ba 3 CuSb 2 09 ^. The simulation temperature, Tf^z, is chosen so as 
to give the best fit to the experimental data and L = 48. (a) Simu¬ 
lated structure factor at T = Tfrz with = tt. (b) Cut through the 
simulated structure factor at qy = ^ and qz = tt (blue dots, shown 
by white dashed line in (a)) compared to x-ray diffraction exper¬ 
iments (red dots). Bragg peaks at qx = ±27r are ignored in the 
simulation, since these are independent of the dumbbell ordering. 


factor is that it can be compared with low temperature 
x-ray diffraction data. This allows the freezing temper¬ 
ature, Tfrz, of the sample to be determined, and then 
simulation at this temperature can be used to shed light 
on the low-temperature structure of the dumbbells in 
the material. One way to determine Tfrz is to consider 
the ratio T = >5(0, 27r/\/3, gz)/5'(27r/3, 27r/\/3, gz), since 
this is sensitive to temperature, as can be seen in Fig. 2. 
The inset to Fig. 2 shows how R increases as a func¬ 
tion of T, eventually saturating in the uncorrelated, high- 
temperature region. 

X-ray diffraction data for Ba 3 CuSb 209 , which is taken 
from Ref. [7], is shown in Fig. 3. The value R ~ 0.4 
is extracted, giving Tfrz/'^nn ~ 0.9, and the simulated 
structure factor at this temperature is superposed on the 
experimental data, showing a good fit. The freezing tem¬ 
perature can be converted into Kelvin by reintroducing 
the dimensionful prefactors in Tcoui [Eq. 1]. The only un¬ 
known is the relative permittivity e^. The dumbbells are 
definitely frozen at T = 300K, the synthesis temperature 
is > lOOOK^, and, for Tfrz to be within these limits, a not 
unreasonable value of ^ 10 is necessary. Furthermore, 
the relatively high value of Tfrz provides a justification 
for ignoring coupling between bilayers. However, the fact 
that diffuse scattering is observed at / = 10 suggests that 
some inter-bilayer correlation is present^. This is left for 
future investigation. 

Once Tfrz has been determined, the dumbbell structure 
at this effective lattice temperature can be studied in 
detail. The density of defect triangles, Utri, on which all 
three dumbbells are orientated in the same direction, is 
shown in Fig. 4. This density is measured relative to a 
ferromagnetic state, in which all dumbbells are orientated 
in the same direction. The density, ntri, increases steadily 
with temperature and, for Tfrz/V^nn = 0.9, is given by 
ntri ~ 0.03. 

Also shown in Fig. 4 is the density of hexagonal plaque- 
ttes, Uhex, measured relative to a long-range honeycomb 
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FIG. 4: Fraction of hexagonal plaquettes, rihexj and defect trian¬ 
gles, TT-tri, as predicted from simulations of Fcoul [Eq. 1]. Hexagonal 
plaquettes have 6 dumbbells of equivalent orientation surrounding 
a dumbbell of the opposite orientation. The fraction of hexagonal 
plaquettes relative to a long range honeycomb lattice {N/3 plaque¬ 
ttes) rapidly saturates with increasing temperature at rihex ~ 0.035 
(blue, upper curve). Defect triangles have three dumbbells with 
the same orientation, and the fraction relative to a ferromagnetic 
state (2N defect triangles) steadily increases with temperature 
(red, lower curve). The black dashed line shows Tj'ipnn = 0.9, which 
is believed to describe the low-temperature dumbbell structure of 
the Ba 3 CuSb 2 09 crystals studied in Ref. [7] (see Fig. 3). 


arrangement of dumbbells {N/3 plaquettes). Hexagonal 
plaquettes are defined as 6 equivalently orientated dumb¬ 
bells surrounding a dumbbell of the opposite orientation. 
The hexagon plaquette density remains low at all temper¬ 
atures, rapidly saturating at only n^ex ~ 0.035, and, for 
Efrz/V^nn = 0.9, is giveu by nhex ~ 0.025. In Ref. [7], the 
presence of diffuse peaks in the x-ray diffraction spectrum 
at = (27r/3, 27r/\/3) [see Fig. 3] was taken as proof of 
a short-range honeycomb arrangement of the dumbbells, 
since this is the wavevector at which Bragg peaks are 
found for a long-range ordered honeycomb arrangement. 
Here we have shown that such a signal arises even in the 
absence of a significant number of hexagonal plaquettes, 
the building blocks of the honeycomb lattice. 

How should the lattice of X ions be described, if not 
by a honeycomb lattice? To answer this a representative 
snapshot of the simulations at Tfrz/'^nn = 0-9 is shown 
in Fig. lb. The lattice can be divided into a set of 
equally orientated clusters - that is clusters of neighbour¬ 
ing dumbbells of the same orientation that are completely 
surrounded by dumbbells of the opposite orientation 
(shown joined by either blue or orange bonds in Fig. lb). 
Superexchange between the electronic degrees of freedom 
associated with the X ions predominantly occurs within 
these equally orientated clusters, as superexchange be¬ 
tween oppositely orientated neighbouring dumbbells is 
expected to be weak^. These superexchange-linked clus¬ 
ters can be seen to have a branching structure, and a 
wide distribution of sizes, n. 

In Fig. 5 we show p(n), the probability that an ar¬ 
bitrary site is part of an n-site cluster. For T = 48 
(N=6912) and for 10 < n < 2000 a good fit to the nu¬ 
merical data is obtained using a power-law probability 
function, p{n) = with C = 0.063 and r = 2.06. 
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FIG. 5: Statistics of superexchange linked Cu^+ clusters at 

Ffrz/t^nn = 0.9, measured by simulation of Fcoul [Eq. 1]. (a) The 
probability, p(n), that a site belongs to a cluster of size n (L = 48). 
For 10 < n < 1000 a power-law distribution p(n) = 0.063/n^-°® 
provides a good fit to the data. At large cluster sizes {n > 1000) 
the finite size of the simulation becomes important, (b) A 6-site 
superexchange linked cluster of Cu^“^ ions, with 5 distinct maximal 
dimer coverings. For geometric reasons, only 4 sites can be covered, 
leaving 2 monomers (orphan sites). 


Finite size effects result in a peak of p{n) at large n and 
the power law also breaks down at n < 10, where stripe¬ 
like correlations between Ising spins suppress the number 
of small clusters. Finite size scaling analysis of the aver¬ 
age size of the largest cluster shows (nmax) oc where 
df = 1.90 is the fractal dimension [see Appendix C for 
more details]. 

These findings match very well to percolation the¬ 
ory for a model of random site filling on the triangu¬ 
lar lattice^^. In this model the percolation threshold is 
at 1/2-filhng, and at this critical point the exponents 
r = 187/91 = 2.055 and df = 91/48 = 1.896 are pre¬ 
dicted. The close agreement between these exponents 
and those found in the simulations suggest that the distri¬ 
bution of sizes of the superexchange linked clusters is at, 
or very close to, the percolation critical point. Thus su¬ 
perexchange linked clusters can be expected at all length- 
scales [see Appendix C]. 

It is common in the Ba 3 XSb 209 family that a sizeable 
fraction of the electronic spins are “orphaned” and in¬ 
teract only weakly with the rest of the system. This is 
observed from a variety of experimental probes and, for 
X=Cu, the percentage of orphan spins has been measured 
in the range 5-16%^“^’^^. Neutron scattering studies pro¬ 
vide evidence that the Cu spins form nearest-neighbour 
singlet bonds at low temperature^, leading us to consider 
covering the lattice in singlet dimers. Maximally cover¬ 
ing the lattice of Cu^+ ions with nearest-neighbour sin¬ 
glet dimers leaves a number of orphan spins, due to the 
geometry of the clusters, and an example of this is shown 
in Fig. 5b. At Tfrz/'^nn = 0.9 the percentage of orphan 
spins calculated in this way is 6%, and, for T > this is 
almost independent of both the simulation temperature 
and the system size [see Appendix D]. 

In this dimer picture, clusters with n = 1 are guar¬ 
anteed to be an orphan spin, and make up about 15% 


of the total orphan spin population. At low tempera¬ 
ture, ESR measures the local environment of the orphan 
spins^, and is therefore biased towards a hexagonal local 
environment. 

It is interesting to speculate about the low tempera¬ 
ture spin-orbital state in Ba 3 CuSb 209 . Theory suggests 
that a nearest-neighbour singlet bond is associated with a 
ferro-orbital alignment between the two sites^^’^^’^^. In 
order for an orbital resonance to occur, it is therefore 
necessary for the system to resonate between different 
singlet coverings of the lattice. The mechanism for this 
resonance can arise directly from the superexchange in¬ 
teraction, or from coupling to the lattice^^’^^’^^. 

For a typical Cu^+ superexchange-linked cluster found 
from solving EqouI [Eq. 1] at T = Tfrz, there are many 
possible maximal dimer coverings, which, for geometri¬ 
cal reasons, leave a number of uncovered monomer sites 
(orphan spins). An oscillation between different dimer 
coverings can equivalently be viewed as a hopping of or¬ 
phan spins around the cluster. Thus resonance between 
different dimer configurations of a cluster not only pro¬ 
vides a mechanism by which orbitals can resonate, but 
also suggests that most orphan spins will be delocalised. 
Since the largest superexchange linked cluster diverges in 
the thermodynamic limit, a resonating state of this type 
can be designated a spin-orbital liquid on the branch lat¬ 
tice. In this picture it is the correlated dumbbell disorder 
that promotes liquid-like behaviour. 

To test this picture we performed exact diagonalisa- 
tion for a spin-orbital Hamiltonian^^ on the 6-site cluster 
shown in Eig. 5 [see Appendix E for details]. A trial 
ground state wavefunction was constructed from the 5 
different dimer coverings of the cluster (shown in Eig. 5), 
and the overlap with the exact wavefunction was 0.98. 

It is feasible to check experimentally the premise of this 
article: that X^+-Sb^+ dumbbells are fiippable at high 
temperature, freeze as temperature is lowered and that 
the low-temperature structure can be understood from 
simulating Ecoui [Eq* 1] at Tfrz- Since Tfrz is controlled by 
the cooling rate, there should be a large difference in the 
dumbbell structure between a crystal slowly cooled from 
the synthesis temperature and one that is quenched, and 
this can be studied by making x-ray diffraction measure¬ 
ments and extracting the ratio R (see Eig. 2). 

In conclusion, we have considered the lattice structure 
of the Ba 3 XSb 209 family, which includes a number of 
proposed spin-liquid materials. By studying a model of 
charged dumbbells on the triangular lattice using Monte 
Carlo simulations, we find a non-trivial lattice structure 
[see Eig. lb], in which superexchange linked clusters of 
X ions form a fractal branching structure. Eocusing in 
particular on X=Cu, which has been proposed as a spin- 
orbital liquid, we show that the obtained lattice structure 
is consistent with x-ray diffraction data. A simple model 
of nearest-neighbour singlet covering of the lattice results 
in a reasonable estimate for the number of orphan spins, 
and gives rise to a scenario in which correlated dumbbell 
disorder promotes a spin-orbital liquid state with non- 
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localised orphan spins. 
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Appendix A: 'ipij{z) interaction matrix 

Here we study in more detail the form of the interaction 
matrix ^pij (z), which appears in Eq. 2 in the main text. In 
particular we consider how this changes with 2 ;, the height 
of the dumbbells [see Eig. la in the main text]. In the 
main text we use 2 : = 0.46a, with a the triangular lattice 
constant, which is taken from the material parameters of 
Ba3CuSb209^. 



FIG. 6: The interaction matrix 'ipij{z)/'iljnn{z) at various values of 
z [see Eq. 2 in the main text]. This matrix takes into account the 
long-range nature of the Coulomb interaction by Ewald summation. 
The size of the system is L = 48 (N = 6912). For 2 ; = 0.01 (red), 
a dipolar interaction is shown as a black dashed line. Deviations 
from the dipolar form at large rij are due to the Ewald summation 
technique, in which interactions between the central cluster and an 
infinite set of image clusters are mapped back onto the interactions 
of the central cluster. 


The energy associated with the interaction of dumb¬ 
bells on sites i and j is given by. 



where the ± depends on the relative orientation and Vij 
is the separation. In the case Vij ^ 2 :, the Taylor ex¬ 
pansion made in the second equality can be truncated at 
first order, and the dumbbell-dumbbell interaction has 
a dipolar form. If 2 ; <C a the dipolar form will be valid 
even for nearest-neighbour interactions, while for ^ a, 
it will only be valid for long-range interactions. 


The interaction matrix 'ipijiz) is calculated from Ecoui 
[Eq. 1 in the main text] by Ewald summation^^. The 
Ewald summation technique allows long-range interac¬ 
tions to be simulated on finite size clusters by tiling the 
infinite plane with a series of images of the central clus¬ 
ter, and taking into account the interactions with these 
image clusters. Therefore ipij{z) includes not only the 
interaction between sites i and j in the original cluster, 
but also interactions of i in the original cluster with j in 
all image tiles. 

The dependence of i^ij{z) on the separation r^j is 
shown in Eig. 6 for various values of z. Eocusing on 
^ = 0.01a, it can be seen that i^ij{z) follows a dipolar 
form at small since the interaction is dominated by 
the central cluster. However, at larger Vij the long range 
interactions become more important with respect to the 
interactions within the central cluster, and the dipolar 
form breaks down. 


Appendix B: Technical details of Monte Carlo 
simulations 

Monte Carlo simulation is used to study Ecoui [Eq. 2 
in the main text]. Here we provide a brief account of the 
technical details. 

A combination of Monte Carlo update methods are 
employed, including single spin flips, parallel tempering 
and a worm algorithm based on mapping the nearest- 
neighbour part of the i^ij{z) interaction [Eq. 2 in the 
main text] onto a loop model on the dual honeycomb 
lattice^^. A typical update step includes 9A attempts 
to flip randomly selected individual spins, 10 calls to the 
worm algorithm and a parallel tempering step. Clusters 
are hexagonal in shape and contain N = lattice sites, 
where L is the length of the hexagon edges. Simulations 
are run using 10^ — 10^ updates, with equal numbers of 
thermalisation and measurement steps. 


Appendix C: Statistics of superexchange linked 
clusters 

Here we consider the superexchange linked clusters of 
X ions, illustrated in Eig. 1 of the main text and stud¬ 
ied in detail in Eig. 5. The relative orientation of the 
neighbouring dumbbell degrees of freedom leads to very 
different superexchange interactions between the X ions. 
If two dumbbells are aligned, the superexchange inter¬ 
action between the associated X ions is comparatively 
strong, and if they are anti-aligned, the interaction is 
weak. This can be seen by comparing the bond angles of 
the X-O-O-X exchange pathways^. In consequence, one 
can think of a set of clusters of superexchange linked X 
ions which, to a first approximation, have no interaction 
with one another. 

Understanding the characteristics of these superex¬ 
change linked clusters is important for understanding the 
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physics of the materials. It is shown in Fig. 4 in the main 
text that the clusters have very few defect-triangles, in 
which all three dumbbells are aligned. Also there are 
very few hexagonal plaquettes, where 6 aligned dumb¬ 
bells surround another dumbbell of opposite orientation. 
This suggests that the clusters are primarily formed of 
1-dimensional chains connected by Y-shape junctions, as 
can be seen in Fig. 1 in the main text. An important 
question to answer concerns the size distribution of these 
clusters: are there many small clusters or is the system 
dominated by a few large clusters? This is equivalent to 
asking whether or not the clusters percolate. 

A useful measure of the cluster size distribution is p{n ), 
the probability that an arbitrary site belongs to a cluster 
of size n. Since every site must belong to a cluster, 

N 

y]p(n) = l, (Cl) 

n=l 

where N is the total number of triangular lattice sites 
in the system. One can also define fN{n), the expected 
frequency of n-site clusters in a system of size Y, as, 

= MU, (C2) 

n 

The average cluster size is given by, 

(^cius) Eti/ivN J:L^Pin)/n ^ ^ 

where (Ydus) is the average number of superexchange 
linked clusters. 

In order to measure p{n) we use Monte Carlo simula¬ 
tion. After each Monte Carlo update we split the sys¬ 
tem into superexchange linked clusters, and determine 
their size. Simulations are carried out for system sizes 
L = 6,12,18,24,36,48, where N = 3L^. The tempera¬ 
ture is set to T = Tfrz = O.Q'^nn, in order to mimic the 
crystals synthesised in Ref. [7]. 

The distribution of cluster sizes for different L is shown 
in Fig. 7. One reason for choosing to measure p{n) is 
that it is independent of Y, and simulations at different 
L collapse onto the same curve. The finite size of the sim¬ 
ulations cuts-off the cluster size at large n, and there is a 
pronounced bump, the position of which is L-dependent. 

It can be seen from Fig. 7 that the power law, 

p{n) = (C4) 

with C = 0.063 and r = 2.06, gives a good fit to the 
Monte Carlo simulations for n > 10. The exponent r 
is known as the Fisher exponent^^. This power law fit 
can also be used to demonstrate that the bump at large 
n is a finite size effect. As an example, for the L = 48 
(Y = 6912) system the fraction of sites which are part of 
a cluster with n > 2000 is given by, 

6912 

p{n) = 0.626, (C5) 

2000 



FIG. 7: Statistics of superexchange linked clusters for different 
system sizes, L, measured by Monte Carlo simulation. The tem¬ 
perature is set to T = 0.9'ipnn- The probability that an arbitrary 
site belongs to an n-site cluster is denoted p(n). For n > 10 this 
follows a power law distribution p{n) = , with r = 2.06, and 

this is in excellent agreement with percolation theory for a model 
of random site filling, which predicts r = 187/91 = 2.055 ... at the 
percolation threshold^^. The power-law behaviour breaks down at 
large n, due to the finite size of the simulations. At n < 10 Ising 
correlations result in a non-power law distribution of cluster sizes. 

which compares well to integrating the power law to in- 
finity, 

POO 

/ Cn^-^ = 0.631. (C6) 

J2OOO 

Thus it can be seen that the bump contains all the clus¬ 
ters that in an infinite system would have n > Y/2, 
but due to the finite size of the simulation are cut-off 
at n < Y/2. 

The average cluster size, (n), can be calculated from 
Eq. C3, and it can be seen in Fig. 8 that it saturates 
with increasing L at (n) = 32.8. More interesting is to 
study the average size of the largest cluster, (rimax), as 
a function of L. As is shown in Fig. 8, this follows a 
power law distribution, (rimax) oc where we measure 
df = 1.90 as the fractal dimension^^. 

The simulation results for r and df can be compared 
to the findings of percolation theory^^. For uncorre¬ 
lated filling of sites in d = 2, theory predicts that at 
the percolation threshold r = 187/91 = 2.055... and 
df = 91/48 = 1.896... in excellent agreement with the 
simulations. For random site filling on the triangular 
lattice, the percolation threshold is at a filling fraction of 

1/2. 

The agreement between the random site filling perco¬ 
lation model and simulations of the Ising system shows 
that correlation between Ising spins is unimportant at 
large distances. This is not surprising, since the temper¬ 
ature is considerably above the Ising critical temperature, 
at which there is a transition to a low-temperature stripe 
ordered state. Since the power law behaviour of p{n) 
breaks down at n < 10, this gives a rough measure of 
the Ising correlation length, ^\s ^ ^/l0. This is a simi- 
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FIG. 8: Average and maximum sizes of superexchange linked clusters, as measured by Monte Carlo simulation, (a) As the system size L is 
increased, the average cluster size saturates at {n) = 32.8. (b) The average size of the largest cluster grows as {rimax} oc , with the fractal 
dimension df = 1.90. This is in good agreement with percolation theory for random site filling, which predicts df = 91/48 = 1.896... at 
the percolation threshold^^. This suggests that the dumbbell model studied here is tuned to the percolation critical point. 
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lar to the correlation length ^is ^ 2 found from fitting a 
Lorentzian function with full width at half maximum of 
1/^is to the Ising structure factor shown in Fig. 2 in the 
main text. 

The power law behaviour of p{n) and (rimax) holds over 
all the simulated system sizes. A minimal conclusion is 
therefore that the percolation correlation length, ^p, is 
large compared to the characteristic lengthscale of the 
largest system size (L=48). Finite size simulations can¬ 
not prove whether the system is exactly tuned to the 
percolation critical point, where ^ oc. What one can 
say is that the ratio of the two dumbbell populations in 
the Ising model is peaked at 1:1, and this peak becomes 
sharper as L is increased. This gives an effective filling 
fraction of 1/2, the critical filling fraction of the random 
site filling percolation model on the triangular lattice^^. 
Thus it appears that the lattice of superexchange linked 
clusters is very naturally tuned to the percolation critical 
point. This result is valid for a large range of tempera¬ 
tures, only breaking down at low-temperature close to 
the Ising critical point, where ^is increases, and Ising cor¬ 
relation effects become important at large lengthscales. 
For temperatures above this, ^\s determines the length- 
scale at which the power law behaviour in p{n) breaks 
down. 


Appendix D: Orphan spins 

Experimentally, it has been shown that Ba 3 CuSb 209 
has an orphan spin population of 5-16%^“^’^^. Here we 
make a theoretical estimate of the number of orphan 
spins that follows from the branch lattice structure, il¬ 
lustrated in Fig. 1 of the main text, and find that this is 
in good agreement with the experimental measurements. 

The branch lattice can be split into superexchange 
linked clusters of dumbbells with the same orientation, as 
described in detail in Section C. To a first approximation, 
superexchange interactions between the X ions are strong 


within the cluster, while the interaction between ions be¬ 
longing to different clusters can be neglected. The simple 
model we propose is that of dimer covering of the clusters. 
In this model dimers correspond to nearest-neighbour 
spin-singlet bonds, with associated ferro-orbital align¬ 
ment. Due to the geometry of the clusters, it is in gen¬ 
eral not possible to cover all sites with dimers, and there 
will remain some monomer sites [see Fig. 5 in the main 
text]. These monomers can be interpreted as orphan 
spins. The idea of nearest-neighbour spin-singlet cov¬ 
ering of the clusters arises from both experiment^’^ and 
theory^^’^^. 

The problem we thus wish to solve is how to maximally 
cover an irregular cluster of sites selected from the trian¬ 
gular lattice with dimers. First we perform Monte Carlo 
simulations of the Ising model given in Eq. 2 of the main 
text. After each Monte Carlo update the lattice is split 
into superexchange linked clusters, and we find a power 
law distribution of cluster sizes, as shown in Eig. C. Eor 
each of these clusters we use a numerical algorithm to 
maximally cover it in dimers, and determine the number 
of monomer sites remaining. 

The algorithm is as follows. Eirst each site of the lat¬ 
tice is assigned a connectivity, which is just the number 
of nearest-neighbour sites on the triangular lattice that 
are part of the same cluster. This connectivity is in the 
range 1 — 6, though in practice connectivities of 4 — 6 are 
relatively uncommon. We select at random one of the 
sites with the lowest connectivity. If this has connectiv¬ 
ity 1 (the usual case) we place a dimer between the chosen 
site and its neighbouring site, delete both sites from the 
cluster and recalculate the connectivitites of the remain¬ 
ing sites. If the smallest connectivity is 2 or greater, a 
dimer is placed on the bond connecting the selected site 
and a randomly chosen second site. These sites are then 
removed from the cluster and the connectivities recalcu¬ 
lated. After the recalculation, it sometimes happens that 
there will be sites with connectivity 0, and these are as¬ 
signed to be monomers. The same process is repeated 







FIG. 9: Fraction of orphan spins, riorph: as measured by Monte Carlo simulation. Superexchange linked clusters are maximally covered 
in dimers, representing nearest-neighbour spin singlets, leaving a number of monomer spins that are considered to be orphans, (a) For 
constant system size, L = 24, the fraction of orphan spins is only weakly temperature dependent, (b) At constant temperature, T = O.9'0nn, 
the fraction of orphan spins is given by riorph = 0.059, independent of the system size, L. 


until all sites have been covered by a dimer or assigned 
as a monomer. 

We do not have a mathematical proof of this algo¬ 
rithm. For small cluster sizes we have tested it against 
the slow but foolproof method of enumerating all dimer 
coverings and finding which one(s) cover the most sites. 
We have also constructed a number of apparently patho¬ 
logical cases, and determined that the algorithm correctly 
calculates the maximal dimer covering. 

The fraction of orphan spins, norph, calculated by this 
method is shown in Fig. 9 as a function of both tem¬ 
perature, T, and of system size L. For fixed system 
size, L = 24, it can be seen that norph varies very lit¬ 
tle with increasing temperature. For fixed temperature, 
T = O.9'0nn5 it can be seen that norph = 0.059 is essen¬ 
tially independent of L. This value is in good agree¬ 
ment with experimental measurements in Ba 3 CuSb 209 
of norph =0.05-0.16^-^’^^. 


where, 

P^=° = I-S,.S„ = (E2) 

are singlet and triplet projection operators and is 
an orbital pseudospin-1/2, with = 1/2 representing 
a orbital and = -1/2 a orbital. The 

vector Uij is different for the 3 possible bond orientations 
(see Fig. 1 in the main text for bond labelling) and given 
by, 

niieA = (0,0,1) 




Appendix E: Spin-orbital state of a 6-site cluster 


In order to test the nature of the spin and orbital state 
of superexchange linked clusters of Cu in Ba 3 CuSb 209 we 
perform exact diagonalisation of a spin-orbital Hamilto¬ 
nian on a 6-site cluster. 

The Hamiltonian we consider was derived in Ref. [12] 
from a 2-orbital Hubbard model, and is given by. 


'Hst 







8 / N 5 ' 

2 16 


} 


UsT [Eq. El] is equivalent to Eq. 8 in Ref. [12] with 
the Hund’s rule coupling set to J = 0 and the ratio of 
hopping parameters set to t/t' = —1/3, which is believed 
to be relevant for Ba 3 CuSb 209 ^^. 



FIG. 10: Trial spin-orbital wavefunction, Ibtrial}? ior Tisj [Eq. El] 
on a 6-site cluster. Blue ellipses represent spin singlet bonds with 
ferro-orbital alignment in the orbitals (and 27r/3 rotations). 

Unpaired sites can be thought of as orphan spins, and are found to 
have a Q 2 type Jahn-Teller distortion (see Fig. 11). The overlap 
with the exact wavefunction, \'iPed): is found to be {'btrialIt^Eo) = 
0.98, showing that the trial wavefunction gives a good description 
of the exact ground state. 
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We perform exact diagonalisation of 1-Lsj [Eq. El] on 
the 6 -site cluster shown in Fig. 10. The idea of choosing 
such a cluster is that it cannot be fully covered with 
dimers, and therefore one can explore the nature of the 
orphan spins (see Section D). 


FIG. 11: Jahn-Teller distortions of the oxygen octahedra associ¬ 
ated with Cg orbitals. A type distortion is expected when the 
orbital state is = ±1/2 (and 27r/3 rotations). A Q 2 type distor¬ 
tion is expected when the orbital is in the = ±1/2 state (and 
27r/3 rotations). 

The exact ground state wavefunction, IV^ed), is found 


to be in the S' = 0 sector. In order to elucidate the physi¬ 
cal nature of this exact wavefunction, we construct a trial 
wavefunction, |'Atrial )5 shown in Fig. 10. In the trial wave- 
function blue ellipses represent nearest-neighbour spin 
singlets with ferro-orbital alignment of the two orbitals in 
either (A-bonds), (B-bonds) or (C- 

bonds). Unpaired sites contain an orphan spin, and it is 
found that the orbital degree of freedom on these sites 
favours a Q 2 type distortion (T^ = ±1/2 and 27r/3 rota¬ 
tions), as opposed to the Qs type distortion (T^ = ± 1/2 
and 27r/3 rotations) favoured by orbitals associated with 
singlet bonds (see Fig. 11 ). We find (V^triailV^ED) = 0.98, 
showing that the trial wavefunction gives a good descrip¬ 
tion of the exact ground state. 

These findings support the picture of the spin-orbital 
state given in the main text. That is, resonating valence 
bonds of nearest-neighbour, ferro-orbital spin singlets 
coexisting with mobile, weakly-coupled orphan spins. 
The finding that orphan spins are associated with a Q 2 
type of Jahn-Teller distortion could be important for un¬ 
derstanding the diffuse x-ray scattering experiments re¬ 
ported in Ref. [9]. 
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In a 3d model of charged dumbbells we expect the ground 


state to retain a stripe pattern within the layers. For 
nearest-neighbour interlayer coupling, the stripe state gives 
the best possible interlayer energy consistent with the ab¬ 
sence of defect triangles in the plane. Other simple 2d 
states, such as the honeycomb, give much worse interlayer 
energies. Thus we expect that the 3d model shows the same 
qualitative features as the 2d model. 


